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Abstract 

The Keller- Segel system has been widely proposed as a model for 
bacterial waves driven by chemotactic processes. Current experiments 
on E. coli have shown precise structure of traveling pulses. We present 
here an alternative mathematical description of traveling pulses at a 
macroscopic scale. This modeling task is complemented with numeri- 
cal simulations in accordance with the experimental observations. Our 
model is derived from an accurate kinetic description of the meso- 
scopic run-and-tumble process performed by bacteria. This model can 
account for recent experimental observations with E. coli. Qualita- 
tive agreements include the asymmetry of the pulse and transition in 
the collective behaviour (clustered motion versus dispersion). In ad- 
dition we can capture quantitatively the main characteristics of the 
pulse such as the speed and the relative size of tails. This work opens 
several experimental and theoretical perspectives. Coefficients at the 
macroscopic level are derived from considerations at the cellular scale. 
For instance the stiffness of the signal integration process turns out to 
have a strong effect on collective motion. Furthermore the bottom-up 
scaling allows to perform preliminary mathematical analysis and write 
efficient numerical schemes. This model is intended as a predictive tool 
for the investigation of bacterial collective motion. 



* University of Edinburgh, School of Mathematics, JCMB, King's Buildings, Edinburgh 
EH9 3JZ, UK. Email: n.bournaveas@ed.ac.uk 

f Institut Curie, UMR CNRS 168 "Physico-Chime-Curie" 11 rue Pierre et 
Marie Curie, F- 75248 Paris cedex 05, France. Emails: axel.buguin@curie.fr, 
j onathan . saragosti@cur ie . f r , pascal . silberzan@curie . f r 

*Ecole Normale Superieure de Lyon, UMR CNRS 5669 "Unite de Mathematiques 
Pures et Appliquees", 46 allee d'ltalie, F-69364 Lyon Cedex 07, France. Email: 
vincent . calvez@ens-lyon . f r 

^Universite Pierre et Marie Curie-Paris 6, UMR CNRS 7598 "Laboratoire Jacques- 
Louis Lions", BC187, 4 place Jussieu, F-75252 Paris cedex 5, and Institut Universitaire 
de France. Email: benoit.perthame@upmc.fr 



1 



1 Introduction 



Since Adler's seminal paper [Tj, several groups have reported the forma- 
tion and the propagation of concentration waves in bacteria suspensions 
[U [32]. Typically, a suspension of swimming bacteria such as E. coli self- 
concentrates in regions where the environment is slightly different such as 
the entry ports of the chamber (more exposed to oxygen) or regions of dif- 
ferent temperatures. After their formation, these high concentration regions 
propagate along the channel, within the suspension. It is commonly admit- 
ted that chemotaxis (motion of cells directed by a chemical signal) is one of 
the key ingredients triggering the formation of these pulses. We refer to [38] 
for a complete review of experimental assays and mathematical approaches 
to model these issues and to [3] for all biological aspects of E. coli. 

Our goal is to derive a macroscopic model for chemotactic pulses based 
on a mesoscopic underlying description (made of kinetic theory adapted 
to the specific run-and-tumble process that bacteria undergo [31] ) . We 
base our approach on recent experimental evidence for traveling pulses (see 
Fig- W)- These traveling pulses possess the following features which we 
are able to recover numerically: constant speed, constant amount of cells, 
short timescale (cell division being negligible), and strong asymmetry in the 
profile. 

We describe as usual the population of bacteria by its density p(t, x) 
(at time t > and position x & R). We restrict our attention to the one- 
dimensional case due to the specific geometry of the channels. The cell den- 
sity follows a drift-diffusion equation, combining brownian diffusion together 
with directed fluxes being the chemotactic contributions. This is coupled 
to reaction-diffusion equations driving the external chemical concentrations. 
In this paper we consider the influence of two chemical species, namely the 
chemoattractant signal S(t,x), and the nutrient N(t,x). Although this is a 
very general framework, it has been shown in close but different conditions 
that glycine can play the role of the chemoattractant [35] . Similarly, glucose 
is presumed to be the nutrient. The exact nature of the chemical species 
has very little influence on our modeling process. In fact there is no need 
to know precisely the mechanisms of signal integration at this stage. The 
model reads as follows: 
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The chemoattractant is assumed to be secreted by the bacteria (at a 
constant rate (3), and is naturally degraded at rate a, whereas the nutrient 
is consumed at rate 7. Both chemical species diffuse with possibly different 
molecular diffusion coefficients. We assume a linear integration of the signal 
at the microscopic scale, resulting in a summation of two independent con- 
tributions for the directed part of the motion expressed by the fluxes us and 
un- We expect that the flux us will contribute to gather the cell density and 
create a pulse. The flux u^ will be responsible for the motion of this pulse 
towards higher nutrient levels. Several systems of this type have been pro- 
posed and the upmost classical is the so-called Keller-Segel equation [22j [28] 
for which fluxes are proportional to the gradient of the chemical: 

u s = x(S)VS, u N = x(N)VN. 

In the absence of nutrient, such systems enhance a positive feedback which 
counteracts dispersion of individuals and may eventually lead to aggregation. 
There is a large amount of literature dealing with this subtle mathematical 
phenomenon (see |34[ [T7] and references therein) . Self- induced chemotaxis 
following the Keller-Segel model has been shown successful for modeling self- 
organization of various cell populations undergoing aggregation (slime mold 
amoebae, bacterial colony,. . . ). In the absence of a chemoattractant S being 
produced internally, this model can be used to describe traveling pulses. 
However it is required that the chemosensitivity function x(iV) is singular 
at N = [22\ . Following the work of Nagai and Ikeda [30j , Horstmann and 
Stevens have constructed a class of such chemotaxis problems which admit 
traveling pulses solutions [19], assuming the consumption of the (nutrient) 
signal together with a singular chemosensitivity. We also refer to [38J for a 
presentation of various contributions to this problem, and to [25J for recent 
developments concerning the stability of traveling waves in some parabolic- 
hyperbolic chemotaxis system. In addition, the contribution of cell division 
to the dynamics of Keller-Segel systems (and specially traveling waves) has 
been considered by many authors (see [241115^ 129] and the references therein). 
However these constraints (including singular chemosensitivity or growth 
terms) seem unreasonable in view of the experimental setting we aim at 
describing. 

An extension of the Keller-Segel model was also proposed in seminal 
paper by Brenner et al. [6] for the self-organization of E. coli. Production 
of the chemoattractant by the bacteria triggers consumption of an exter- 
nal field (namely the succinate). Their objective is to accurately describe 
aggregation of bacteria along rings or spots, as observed in earlier exper- 
iments by Budrene and Berg {B\ that were performed over the surface of 
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gels. One phase of the analysis consists in resolving a traveling ansatz for 
the motion of those bacterial rings. However the simple scenario they first 
adopt cannot resolve the propagation of traveling pulses. The authors give 
subsequently two possible directions of modeling: either observed traveling 
rings are transient, or they result from a switch in metabolism far behind 
the front. The experimental setting we are based on is quite different from 
Budrene and Berg's experiments (in particular regarding the dynamics): for 
the experiments discussed in the present paper, the bacteria swim in a liquid 
medium and not on agar plates. Therefore we will not follow [6]. 

On the other hand Salman et al. [35] consider a very similar experimental 
setting. However the model they introduce to account for their observations 
is not expected to exhibit pulse waves (although the mathematical analysis 
would be more complex in its entire form than in [H]). Actually Fig. 5 
in [35] is not compatible with a traveling pulse ansatz (because the pulse 
amplitude is increasing for the time of numerical experiments). 

Traveling bands have also been reported in other cell species, and espe- 
cially the slime mold Dictyostelium discoideum [39] . Notice that the original 
model by Keller and Segel [22] was indeed motivated by the observation of 
traveling pulses in Dictyostelium population under starvation. This question 
has been developped more recently by Hofer et al. [18] using the Keller- 
Segel model, as well as Dolak and Schmeiser [11] and Erban and Othmer 
[13] using kinetic equations for chemotaxis. According to these models, the 
propagating pulse waves of chemoattractant (namely cAMP) are sustained 
by an excitable medium. The cells respond chemotactically to these waves 
by moving up to the gradient of cAMP. Great efforts have been successfully 
performed to resolve the " back-of-the-wave paradox": the polarized cells are 
supposed not to turn back when the front passed away (this would result 
in a net motion outwards the pulsatile centers of cAMP). Although we are 
also focusing on the description of pulse waves, the medium is not expected 
to be excitable and the bacteria are not polarized. Nevertheless, we will re- 
tain from these approaches the kinetic description originally due to Alt and 
coauthors [21 [31] . This mathematical framework is well-suited for describing 
bacterial motion following a microscopic run-and-tumble process. 

A new class of models for the collective motion of bacteria has emerged 
recently. It differs significantly from the classical Keller-Segel model. Rather 
than following intuitive rules (or first order approximations), the chemotac- 
tic fluxes are derived analytically from a mesoscopic description of the run- 
and-tumble dynamics at the individual level and possibly internal molecular 
pathways, see pH QUI H QU [33] • The scaling limit which links the 
macroscopic flux us (or similarly um) to the kinetic description is now well 
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Figure 1: Experimental evidence for pulses of Escherichia coli traveling across a 
channel. The propagation speed is constant and the shape of the pulse front is 
remarkably well conserved. Observe that the profile is clearly asymmetric, being 
stiffer at the back of the front (see also Fig. [2]). Cell division may not play a crucial 
role regarding the short time scale. 



understood since the pioneering works [21 [31]. Here we propose to follow the 
analysis in [11], which is based on the temporal response of bacteria [71 136]. 
denoted by (fig in Appendix ()A.8p , Namely we write these fluxes as: 

us = Xs.U(-ed t S,\VS\) — , (1.2) 

where e is a (small) parameter issued from the microscopic description of 
motion. Namely e is the ratio between the pulse speed and the speed of 
individual cells (they differ by one order of magnitude at least according 
to experimental measurements). The function contains the microscopic 
features that stem from the precise response of a bacterium to a change 
in the environment (see the Appendix ()A.8|) ) . It mainly results from the 
so-called 'response function' at the kinetic level that describe how a single 
bacterium responds to a change in the concentration of the chemoattractant 
S in its surrounding environment. We give below analytical and numerical 
evidence that traveling pulses exist following such a modeling framework. 
We also investigate the characteristic features of those traveling pulses at 
the light of experimental observations. 

The experiments presented in the present paper will be described in 
more details in a subsequent paper. Briefly, in a setup placed under a low 
magnification fluorescence microscope maintained at 30 °C, we fill polymer 
microchannels - section: ca. 500^m x lOOfim, length ca. 1cm - with a sus- 
pension of fluorescent E. coli bacteria - strain RP437 considered wild-type 
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for motility and chemotaxis, transformed with a pZElR-gfp plasmid allow- 
ing quantitative measurement of bacteria concentration inside the channel. 
We concentrate the cells at the extremity of the channel and monitor the 
progression of the subsequent concentration wave along the channel. In par- 
ticular, we dynamically extract the shape of the front and its velocity (see 
Fig. I]). 

Coupling the model (jl.lB with the formula (|1.2p results in a parabolic 
type partial differential equation for the bacterial density p, such as in the 
Keller-Segel system. It significantly differs from it however, as it derives 
in our case from a kinetic description of motion. Especially the flux us 
is uniformly bounded, whereas the chemotactic flux in Keller-Segel model 
generally becomes unbounded when aggregative instability occurs, which is 
a strong obstacle to the existence of traveling pulses. 



2 Traveling pulses under competing fluxes 

2.1 Stiff response function: pulse wave analytical solutions 

It is usually impossible to compute explicitely traveling pulse solutions for 
general systems such as (jl.ip . To obtain qualitative properties is also a 
difficult problem: we refer to (30J US EH] for examples of rigorous results. 
Here, we are able to handle analytical computations in the limiting case of 
a stiff signal response function, when the fluxes are given by the expression 
(see the Appendix (|A,8|>): 

2^ 

US = XS\1- l^l I sign(d x S), (2.1) 




u N = XN 1 - ^7 si g n (W • ( 2 - 2 ) 



In other words, a specific expression for Ja, in (II. 2p is considered in this 
section. It eventually reduces to us = xssign (d x S) as e — ► 0. 

We seek traveling pulses, in other words particular solutions of the form 
p{t,x) = p(x — at), S(t,x) = S(x — at), N(t,x) = N(x — at) where a 
denotes the speed of the wave. This reduces (jl.ip to a new system with a 
single variable z = x — at: 

{ -ap'(z) = D p p"(z)-(p(z)u s (z)+p(z)u N (z))' , 

-aS'(z) = D s S"{z)-aS(z)+l3p(z), (2.3) 
-aN'(z) = D N N"{z)-~fp{z)N{z). 
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We prescribe the following conditions at infinity: 



p(±oo) = 



S(±oo) = 



N(±oo) = N± . 



(2.4) 



We impose a > without loss of generality. This means that the fresh 
nutrient is located on the right side, and thus we look for an increasing 
nutrient concentration N'(z) > 0. We expect that the chemoattractant 
profile exhibits a maximum coinciding with the cell density peak (say at 
z = 0), and we look for a solution where S'(z) changes sign only once at 
z = 0. Then, the fluxes (I2,ip - (|2.2p express under the traveling wave ansatz 
as: 



The flux us takes two values (with a jump at z = 0), whereas the flux um 
is constant. Therefore the cell density is a combination of two exponential 
distributions: 



V up 

This combination of two exponentials perfectly match with the numerical 
simulations (Fig. [5]). 

2.2 Formula for the traveling speed 

To close the analysis it remains to derive the two unknowns: the maximum 
cell density po an d the speed a, given the mass and the constraint that d z S 
vanishes at z = (because S(z) reaches a maximum). 

On the one hand, the total mass of bacteria is given by M = po(l/\~ + 
1/|A + |). On the other hand the chemotactic field is given by S(z) = (K * 
/3p)(z), where the fudamental solution of the equation for S(z) is 




Integrating once the cell density equation we obtain 



d pP'{z) = P{ z ) (us(z) + u N (z) - a) . 
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Figure 2: (Top) Experimental evidence of a traveling pulse: time snapshots of the 
full experiments described in Fig. [TJ The density profile is clearly asymmetric and 
preserved along the time course of the experiment. Main contribution to growth 
takes place at the back of the pulse. This suggests that nutrients are not totally 
consumed by the pulse. The number of bacteria in the pulse is approximately 
constant during the pulse course. (Bottom) A generic density profile obtained 
with the model (jl.lj) (see also Fig. [3|). The bacteria (plain line) are attracted by 
a nutrient (dashed line) which is consumed. They secrete in addition their own 
chemical signal (dash-dotted line) which is concentrated at the peak location. 
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To match the transition in monotonicity condition, the chemical signal 
should satisfy S'(0) = 0, that is (K' * /3p)(0) = 0, which leads to 



po 



= oi 



f 

/ a\{a 2 + a 3 ) exp(a2-z + a 3 z) exp(A z) dz 

J — oo 

+Po / ai(-a 2 + a 3 ) ex.p(-a 2 z + a 3 z) exp(\ + z) dz 
Jo 



a 2 + 03 — a 2 + 03 



a 2 + 03 + —a>2 + a 3 + \ + 



This leads to the following equation that we shall invert to obtain the front 
speed: 

A~ a 2 + a 3 



\X+\ a 2 — a 3 

\2 X 



a + Us + Xn) [1 ~ (ea) ) ^AD s a + a 2 + a 



From this relation we infer: 



(2.5) 



a a 

= X.S 

1 - {eaf ^ AD sot + a 2 



~ ~ — ; — 72 = xs /m „ == • (2.6) 



We deduce from monotonicity arguments that there is a unique possible 
traveling speed a £ (0,e _1 ). 

According to f|2.6|) the expected pulse speed does not depend upon the 
total number of cells when the response function is stiff. This can be related 
to a recent work by Mittal et al. [27] where the authors observe experimen- 
tally such a fact in a different context (see Section T2.3I below) . In the case of 
a smooth tumbling kernel in (|A.8jl . our model would predict a dependency 
of the speed upon the quantity of cells. But this analysis suggests that the 
number of cells is presumably not a sensitive biophysical parameter. 

Observe that the speed does not depend on the bacterial diffusion co- 
efficient either. Therefore we expect to get the same formula if we follow 
the hyperbolic approach of [11] in order to derive a macroscopic model. 
Indeed the main difference at the macroscopic level lies in the diffusion co- 
efficient which is very small in the hyperbolic system. Nevertheless, the 
density distribution would be very different, being much more confined in 
the hyperbolic system. Furthermore, scaling back the system to its original 
variables, we would obtain a pulse speed being comparable to the individual 
speed of bacteria (see Appendix). This is clearly not the case. 
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2.3 Cluster formation 



Stiff response function. 

Mittal et al. have presented remarkable experiments where bacteria E. coli 
self-organize in coherent aggregated structures due to chemotaxis |27j . The 
cluster diameters are shown essentially not to depend on the quantity of 
cells being trapped. This experimental observation can be recovered from 
direct numerical simulations of random walks [21 j . 

We can recover this feature in our context using a model similar to (jl.ip 
derived from a kinetic description. Following Section 12.11 we compute the 
solutions of (|2.3|) in the absence of nutrient (assuming again a stiff response 
function). Observe that stationary solutions correspond here to zero-speed 
traveling pulses. The problem is reduced to finding solutions of the following 
system: 



We assume again that sign (S'(x)) = — sign(x). This simply leads to, 



This is compatible with the postulate that S(x) changes sign only once, at 
x = (the source (3p(x) being even). The typical size of the clusters is of 
the order A -1 , which does not depend on the total number of cells. This 
is in good agreement with experiments exhibited in [27] . The fact that we 
can recover them from numerical simulations indicates that these stationary 
states are expected to be stable. 

General response function. 

Cluster formation provides a good framework for investigating the situation 
where we relax the stiffness assumption of the response function cpg. Below 
(fis is characterized by the stiffness parameter 5 through <^(0) = —3/5 (see 
Appendix). Consider the caricatural model (in nondimensional form): 




(2.7) 



p(x) = po exp(— A|x|) , where A = — — . 

Up 



< 




p-d x (pu s ) , u s 



2 



1 



/ 



»e(-i,i)) 



v(frs(edtS + vd x S) dv 



d xx S + aS = p . 



(2.8) 
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We rewrite a = l~ 2 , where I denotes the range of action of the chemical 
signal. We investigate the linear stability of the constant stationary state 
(p, S) = ((p), a^ 1 {p)) where (p) is the meanvalue over the domain [0, L]. We 
introduce the deviation to the stationary state: n = p — ~p, c = S — S. Then 
the linearized system writes close to (~p,S): 

d tn = d 2 xx n - d x ((p)u s ) , u s = -- / v 2 (p' s (0)d x cdv = , 

—d 2 x c + ac = n . 

(2.9) 

The associated eigenvalue problem reduces to the following dispersion rela- 
tion for £ = 27rk/L, 

m = -e + M e 



5La + i 2 ' 

Due to the conservation of mass, we shall only consider k > 1. The eigen- 
value becomes positive if there exists k such that 

M 2 ,,^^£, 2271 

->a + e, or equivalents — > ^ + k - . 

Therefore the constant solution is linearly stable if the following condition 
is fulfilled: 

ML 1 /L x2 



2vr5 2vr V I 

The picture is not complete as we have not investigated the stability of 
the non-trivial steady-state. However this indicates that the mass and the 
stiffness parameter 5 play important roles regarding cluster formation. 



3 Numerical insights 

We complete the theoretical analysis with some numerical simulations ex- 
hibiting propagation of pulses (or not) in regimes where analytical solutions 
are not available (see Fig. [3]). 

We opt for the following initial conditions in our numerical experiments: 
a decreasing exponential function centered on the left side of the channel for 
the cell density, no chemical signal, and a constant level of nutriment Nq. 

Parameters are issued from literature (see [35] and references therein) 
and from the mesoscopic derivation of the model (see Appendix 'Parameters 
and scales'). Adimensionalizing time and space yields the time and space 
scales being respectively 10s and 200pm, such that the total duration of the 
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Figure 3: (Top) Experimental results under abundant nutrient conditions: M9 
minimal medium supplemented with 4% glucose and 1% casamino acids (both ten 
times more concentrated than in the case of Fig. [5}. (Bottom) Numerical simula- 
tions of system ([l.ip in the case of unlimited nutrient, and a stiff response function 
(j>s- We observe the propagation of a traveling pulse with constant speed and asym- 
metric profile. The value of the speed and the features of the profile (combination of 
two exponential tails) perfectly match the traveling wave ansatz analysis of Section 
[2j We use a semi-implicit upwind finite-difference scheme performing a half-point 
discretization (in space) of the time derivative of chemical species when comput- 
ing the approximate fluxes. Specific parameters are: 5 = 10 -3 and N = 10 (see 
Section [3] for the other parameters). 
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Figure 4: (Top) In this experiment, bacteria are cultivated at a concentration of 
5. 10 8 cells/ml in the same rich medium as in Fig. 3. After, they are resuspended 
in LB nutrient to an OD600 of 3. 10 8 cells/ml. We interpret the absence of pulse 
propagation as following. Bacteria are adapted to a rich environmnent before re- 
suspension. Thus they are not able to sense small chemical fluctuations necessary 
for clustering to occur when evolving in a relatively poor medium. (Bottom) In- 
fluence of the internal processes stiffness. When the individual response function 
4>s is not stiff, the effect of dispersion is too strong and no pulse wave propagates, 
as opposed to Fig. [3j Specific parameters are: 5 = 10 _1 and No = 10 (see Section 
[3] for the other parameters). In mathematical models of bacterial chemotaxis, it is 
commonly accepted that adaptation of cells to large chemoattractant changes acts 
through the measurement of relative time variations: S~ 1 DS/Dt (see Appendix). 
In our context, this is to say that the stiffness parameter 6 is proportional to the 
chemical level S. Hence after having dramatically changed the environment and 
before bacteria adapt themselves, we can consider that the response function (f)$ is 
not stiff. , 
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Figure 5: At low level of nutrient the cell population splits into two subpopula- 
tions. A fraction remains trapped at the boundary (as a stationary profile) and a 
fraction travels accross the channel with constant speed (see also Fig. Specific 
parameters arc: S = 1CP 3 and Nq = 1 (see Section [3] for the other parameters). 
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Figure 6: (Left) Respective profiles for the bacterial density (plain line), the chem- 
ical concentration (dash-dotted line) and the nutrient concentration (dashed line) 
in the limited nutrient/stiff response function configuration described in Fig. [5] 
(5 = 1CP 3 and Nq = 1). (Right) The net chemoattractive flux us + is plotted. 
It clearly shows a transition between the traveling wave (positive flux) and the 
equilibrium (negative flux) located at the left boundary. 
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computation is approx. lh, and the length of the computational channel is 
4cm. We fix the following parameters, D p = 1, Xs = Xn = 1> for the cell 
density equation; D$ = 2, a = 0.05 for the chemical concentration [35]. We 
assume in addition = for the sake of simplicity: the nutriment is not 
required to act as a non-local signal (and we expect this coefficient to have 
very little influence on the dynamics). By adimensionalizing further the sys- 
tem, we may choose M = 1, 7 = /3 = 1 without loss of generality. Finally the 
signal response function <j>g is chosen as follows: (f>$(Y) = — 2/ir arctan(Y/<5), 
with a stiffness factor 5. We also keep the memory of the drift-diffusion limit 
performed in the Appendix, by setting e = 0.1. The only free parameters 
subject to variations are Nq and 5. 

We can draw the following main conclusions from numerical experiments. 

Influence of the stiffness of the internal response function (Fig. [4]). 

When the stiffness assumption for the internal response function is relaxed, 
no pulse propagation is observed numerically. Deriving the exact conditions 
that guarantee the propagation of a traveling pulse seems to be a challeng- 
ing task. However we give below some heuristics for the particular choice 
4>s(Y) = 5 being a stiffness factor. 

Although the chemotactic equation of (jl.ip is significantly different from 
the Keller-Segel model, they coincide as far as the stability of the homoge- 
neous (unclustered) configuration is under question. We learn from Section 
12.31 that the stiffness parameter 5 plays an important role in the stability 
of the homogeneous solution. It is well known that the Keller-Segel system 
is subject to a bifurcation phenomenon due to its quadratic, non-local non- 
linearity. This is well understood in two dimensions of space for instance 
[341 117] . If some nondimensional parameter is small enough, diffusion dom- 
inates and no self-organization arises in the system. On the contrary, self- 
attraction between cells overcomes diffusion when this parameter is large, 
and yields the formation of a singularity (i.e. aggregation point) [34J. 

Clearly the same kind of mechanism acts here (see Fig. [3] as opposed to 
Fig. [4]). However there is no mathematical argumentation to sustain those 
numerical and intuitive evidence yet. 

Limited versus non limited nutrient (Fig. [5]). 

When the nutrient is limited in the experimental device (and conditions for 
a pulse to travel are fulfilled) then only part of the bacterial population 
leaves the initial bump. The solution seems to be the superposition of a 
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traveling pulse and a stationary state (such as in Section 12,31 in the absence 
of nutrient). Solitary modes with smaller amplitudes may appear at the 
back of the leading one (not shown). To predict which fraction of mass 
starts traveling turns out to be a difficult question. 

4 Conclusions and perspectives 

We present in this article a simple mathematical description for the collective 
motion of bacterial pulses with constant speed and asymmetric profile in a 
channel. The nature of this model significantly differs from the classical 
Keller-Segel system although it belongs to the same class of drift-diffusion 
equations. Our model is formally derived from a mesoscopic description of 
the bacterial density, which allows for a more accurate expression of the 
cell flux. We are able to compute analytically the speed of the pulse and 
its profile in the limit of a stiff response function <f>g. The theoretical pulse 
speed has some striking features: it does not depend on the total number of 
bacteria, neither on the bacterial diffusion coefficient. This can be related 
to experimental evidence by Mittal et al. [27] where bacteria self-organize 
into size-independent clusters. Our approach can be summarized as follows: 
a nutrient is added to pull chemotactic clusters of cells. This creates an 
imbalance in the fluxes which induces the asymmetry of the traveling profile. 

The next step would be to work at the kinetic level. Much has to be done 
for the design of efficient kinetic schemes for the collective motion of cells 
subject to chemotactic interactions. It would also be feasible to point out 
the dependency of the tumbling operator upon some internal variable {e.g. 
the cytoplasmic concentration of protein CheY). This approach carries out 
the coupling between an internal protein network and the external chemoat- 
tractant signals [371 [2Tj . Kinetic models are also relevant for describing this 
microscopic mechanism [121 [5] (the network is basically transported along 
the cells' trajectories). However the increase in complexity forces to reduce 
the size of the network, and to use rather caricatural systems mimicking 
high sensitivity to small temporal variations (excitation) and adaptation to 
constant levels of the chemoattractant. 

Assuming independent integration of the chemical signals constitutes a 
strong hypothesis of our model. There exist two main membranous receptors 
triggering chemotaxis, namely Tar and Tsr. As the signals which act in 
the present experiments are not perfectly determined, we have considered 
the simplest configuration. To further analyse the interaction between the 
external signals, one should include more in-depth biological description of 
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the competition for a single class of receptor [23J . 



A Appendix: Kinetic models for chemotaxis and 
their drift-diffusion limit 

Kinetic framework. 

The classical theory of drift-diffusion limit for kinetic modeling of bacterial 
chemotaxis is a way to compute the macroscopic fluxes ug, ujy in (jl.lj) 
|16j . Because we assume a linear integration of the different signals for each 
individual, we restrict ourselves to the action of a single chemical species S. 

The kinetic framework is as follows. A population of bacteria can be 
described at the mesoscopic scale by its local density f(t,x, v) of cells located 
at the position x and with velocity v at time t. The kinetic equation proposed 
in the pioneering works of Alt, Dunbar and Othmer [21 [31] combines free 
runs at speed v, and tumbling events changing velocity from v' (anterior) to 
v (posterior), resp. according to the Boltzman type equation: 



The velocity space V is bounded and symmetric, usually V = B(0,c) or 
V = S(0, c) (bacteria having presumably constant speed). As we deal with 
the idealization of a two-dimensional phenomenon in one dimension of space, 
we shall perform our computations for V = [— c, c], but the results contained 
in this paper do not depend on this particular choice. Kinetic models of 
chemotaxis have been studied recently in [201 [3[ [5]. 

The turning kernel T describes the frequency of changing trajectories, 
from v' to v. It expresses the way external chemicals may influence cell tra- 
jectories. A single bacterium is able to sense time variations of a chemical 
along its trajectory (through a time convolution whose kernel is well de- 
scribed since the experiments performed by Segall et al. [36]). For the sake 
of simplicity we neglect any memory effect, and we assume that a cell is able 
of sensing the variation of the chemical concentration along its trajectory. 
Following [TT], this is to say that T is given by the expression 




(A.l) 



T[S](v' - v) = iP — = ip {8 t S + v' ■ V X S) . 




(A.2) 
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Figure 7: Tumbling frequencies (at the mesoscopic scale) obtained from the nu- 
merical experiment described in Fig. [SJ the tumbling probability is higher when 
moving to the left (upper dashed line) at the back of the pulse, whereas the tum- 
bling probability when moving to the right is lower (upper plain line), resulting in 
a net flux towards the right, as the pulse travels (see Fig. [5]). Notice that these two 
curves are not symmetric w.r.t. to the basal rate 1, but the symmetry defect is of 
lower order (1CP 3 ). The peak location is also shown for the sake of completeness 
(lower plain line). 



The signal integration function ip is non-negative and decreasing, expressing 
that cells are less likely to tumble (thus perform longer runs) when the 
external chemical signal increases (see Fig. [7] for such a tumbling kernel 
in the context of the present application). It is expected to have a stiff 
transition at 0, when the directional time derivative of the signal changes 
sign [36], E3 EI]- Our study in Section [3] boils down to the influence of 
the stiffness, by introducing a one parameter family of functions ipsiY) = 



Parameters and scales. 

The main parameters of the model are the total number of bacteria M which 
is conserved, the maximum speed of a single bacterium c = max{|v|; v G V}, 
and the mean turning frequency Ao = ipoc d (where d denotes the dimension 
of space according to our discussion above). The main unknown is the speed 
of the traveling pulse, denoted by a. We rescale the kinetic model (jA.ip into 
a nondimensional form as follows: 

t = ti , x = xx , v = v c, V = cV , T = TipQ . 
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We aim at describing traveling pulses in the regime x = at. Experimental 
evidence show that the bulk velocity a is much lower than the speed of a 
single bacterium c. This motivates to introduce the ratio e = a/c. According 
to experimental measurements, we have e ~ 0.1. The kinetic equation writes: 

ed-J + v ■ Vj = ^ ( / i> Ld~ t S + v' ■ VzS) f(t, x, v') dv' 

-\V$(ed { S + v-V z s)f(t,x,v)} , (A.3) 

where ip{z) = i(j(cz/x). Following the experimental setting (see Introduc- 
tion, Fig. [T] and Fig. [2D and the biological knowledge [3], we choose the 
scales x « 200^m, Ao ~ Is , and c = 20/xm.s _1 . Hence \qx/c rs 10. 
Therefore we rewrite this ratio as: 

\qX fJ, 

c e 

where the nondimensional coefficient [i is of order 1. 
Drift-diffusion limit of kinetic models. 

To perform a drift-diffusion limit when e — > (see [HI [TOJ |9l 151] , and 
[111 1 14] for other scaling limits, e.g. hyperbolic), we shall assume that the 
variations of ij) around its meanvalue ipo are of amplitude e at most. It writes 
in the nondimensional version as follows: ipiY) = 1 + etpsiY). Hence the 
chemotactic contribution is a perturbation of order e of a unbiased process 
which is constant in our case because the turning kernel does not depend 
on the posterior velocity and the first order contribution is required to be 
symmetric with respect to (v',v). This hypothesis is in agreement with 
early biological measurements. It is also relevant from the mathematical 
viewpoint as we are looking for a traveling pulse regime where the speed 
of the expected pulse is much slower than the speed of a single individual. 
This argues in favour of a parabolic scaling as performed in this Appendix. 

The rest of this Appendix is devoted to the derivation of the Keller-Segel 
type model in one dimension of space: 

d t p + d x (-D p d x p + pu s ) =0. (A.4) 

Dislike the classical Keller-Segel model (used for instance by Salman et al. 
[35]). singularities cannot form (excessively populated aggregates) with the 
chemotactic flux us given in (|A.8|) below. This is because the latter re- 
mains uniformly bounded (see also Mittal et al. |27] where clusters emerge 
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which are plateaux and thus not as singular as described for KS system in 
a mathematical sense). 

Sketch of parabolic derivation. We start from the nondimensional kinetic 
equation (|A.3|1 : 

ed t f + vV x f = ^\ [ {l + e<t> s [S](v'))f(t,x,v')dv' 
e {.Jv'ev 

-\V\(l + eUS)(v))f(t,x,v)} , 

which reads as follows, 

ed t f + v-V x f = ^ (p(t, x) - \V\f(t, x, v)) 

+ I <t> s [S](v')f(t,x,v')dv' -\V\MS](v)f(t,x,v)) . (A.5) 

Therefore the dominant contribution in the tumbling operator is a relaxation 
towards a uniform distribution in velocity at each position: f(t,x,v) = 
p(t,x)F(v) as e — > 0, where F(v) = \V\ lr ve y\ ■ Notice that more involved 
velocity profiles can be handled [10j[33], but this is irrelevant in our setting 
as the tumbling frequency does not depend on the posterior velocity v. 

The space density p(t, x) remains to be determined. For this purpose 
we first integrate with respect to velocity v and we obtain the equation of 
motion for the local density p(t,x) = j v&v f(t,x,v)dv: 

dtP + ^-j = 0, i = e -1 / vf(t,x,v)dv. 

Jv&V 

To determine the bacterial flow j we integrate (|A.5p against v: 

tdtl / vf(t,x,v)dv) + V x - I / v ®vf{t,x,v)dv\ 
\Jvav J \Jvev J 



= -p\V\j-p\V\ / vMS](v)f(t,x,v)dv. 
Jvev 

We obtain formally, as e — > 0: 

J = -V, [j-\ 2 dv) -p(t,x)±j^vMSKv)dv. 

(A.6) 
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Finally, the drift-diffusion limit equation (|A.4|) reads in one dimension of 
space: 

%P =(7- ! \v\ 2 dv) d 2 xxP + d x (p [ v<p s (ed t S + vd x S) ^ j . 

Jve[-i,i] J \ Jve[-i,i] 2 J 

(A.7) 
□ 

To sum up, we have derived a macroscopic drift-diffusion equation, where 
the bacterial diffusion coefficient and the chemotactic flux are given by: 

D P = -^[ \v\ 2 dv, u s = - [ v<j) 5 {tdtS + vd x S) (A.8) 
4 ^ Jve[-i,i] Jve[-l,i] 2 

In the limiting case where the internal response function (f>s is bivaluated: 
4>s(Y) = ^ol{y<o} — <^ol{y>o} ) the flux rewrites simply as: 



<n / , / d t S^ " 



For the sake of comparison, we highlight the corresponding expressions 
which have been obtained by Dolak and Schmeiser. In [11] authors perform 
a hyperbolic scaling limit leading to the following chemotactic equation for 
the density of bacteria: 

d tP + V- (-eDVp + pUs) =0, 
where D is an anisotropic diffusion tensor and the chemotactic flux is given 



Aj veV if>(dtS + vi7S\) y|VS|' 

for some renormalizing factor A. The two approaches do not differ that 
much at first glance (especially when ip is bivaluated). Notice however that 
the "small" e parameter does not appear at the same location: in front of 
the diffusion coefficient in the hyperbolic limit and inside the chemotactic 
flux in the parabolic limit. 
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